// THIS FILE FINDS THE LAMBDAS (THE MC OF THE MARGINAL PLANT)
// CALCULATES THE AVERAGE MB (TABLE 1)

set more off
cd "${PATH}/data"


use "entsoe.dta", clear

drop eventtime // use time_utc for remaining analysis


// KEEP OBSERVATIONS FOR 2015 and 2016 ONLY	
drop if yofd(daily) > 2016


// Ajustment for fossiloil: TenneT has one plant running at 1 MWH capacity always from March 2016 onwards 
* Do only consider as marginal source if prodcution >1
replace fossiloil = 0 if fossiloil == 1

******************************************
// Generate table with total output per technology and year
preserve
local fuel_types "windoffshore windonshore solar otherrenewable marine hydrowaterreservoir hydrorunofriverandpoundage hydropumpedstorage geothermal waste biomass nuclear fossilbrowncoallignite fossilcoalderivedgas fossilhardcoal fossilgas fossilpeat fossiloilshale fossiloil"

// production measured in MW - construct MWh by aggregating data at the hourly level
collapse (mean) `fuel_types' year (first) TSO , by(daily hour) // use averages 

encode TSO, gen(TSO_num) 
collapse (sum) `fuel_types' , by(TSO year) 
// Total production in MWh

egen fuel_hardcoal = rowtotal(fossilcoalderivedgas fossilhardcoal)
gen fuel_lignite = fossilbrowncoallignite
gen fuel_gas = fossilgas
gen fuel_oil = fossiloil

keep fuel_* TSO year

reshape long fuel_ , i(year TSO) j(type) string
ren fuel_ productionMWh
ren type fuel
save "Production_TSO_annual_fossil.dta", replace
restore
******************************************


// SORT GENERATION SOURCES WITH INCREASING COST
order windoffshore windonshore solar otherrenewable marine hydrowaterreservoir hydrorunofriverandpoundage hydropumpedstorage geothermal waste biomass nuclear fossilbrowncoallignite fossilcoalderivedgas fossilhardcoal fossilgas fossilpeat fossiloilshale fossiloil 


// FIND MARGINAL SOURCE 
local ele_type "windoffshore windonshore solar otherrenewable marine hydrowaterreservoir hydrorunofriverandpoundage hydropumpedstorage geothermal waste biomass nuclear fossilbrowncoallignite fossilcoalderivedgas fossilhardcoal fossilgas fossilpeat fossiloilshale fossiloil" 
local types: word count `ele_type'

forvalues n = 1/`types' {
	local source: word `n' of `ele_type'  
	gen marg_`n' = "`source'" if !missing(`source') & `source'!=0
} 

gen marginal = ""
foreach n of numlist 19(-1)1 {
	replace marginal = marg_`n' if !missing(marg_`n') & missing(marginal)
}

// 50 HZ, special case: Adjustment in fossiloil  //
// REPLACE MARGINAL TECHNOLOGY WITH THE SECOND MOST EXPENSIVE IN 50HZ as OIL production linked to refinery (always running and not responding to prices)
replace marginal = "" if TSO == "50Hertz"
foreach n of numlist 18(-1)1 {
	replace marginal = marg_`n' if !missing(marg_`n') & missing(marginal)  & TSO =="50Hertz"
}
*

drop marg_*


// MARGINAL COSTS ($/MWH)
gen marginal_cost = 0 
replace marginal_cost = 5.1 if marginal == "nuclear" // source: ENTSO-e TYNDP 2018 modeling data as in Golub et al. 
replace marginal_cost = 6 if marginal == "biomass" // allow biomass to be marginal: assign positive value higher than nuclear in line with industry reports
replace marginal_cost = 11.3 if marginal == "fossilbrowncoallignite"  // Brown Coal Report: est. cost of extraction in Germany
replace marginal_cost = cost_MWH_coal if marginal == "fossilhardcoal" // Based on international markets
replace marginal_cost = cost_MWH_ng if marginal == "fossilgas" // Based on international markets
replace marginal_cost = cost_MWH_crude_oil if marginal == "fossiloil" // Based on international markets


// MARGINAL EMISSIONS (t/MWH)
gen marginal_emissions = 0 
replace marginal_emissions = emission_MWH_coal if marginal == "fossilhardcoal"
replace marginal_emissions = 1.148 if marginal == "fossilbrowncoallignite" 
replace marginal_emissions = emission_MWH_ng if marginal == "fossilgas" 
replace marginal_emissions = emission_MWH_crude_oil if marginal == "fossiloil"


**********
// AS DO NOT CONSIDER IMP/EXP, RECREATE TOTAL LOAD AS "ALL ELECTRICITY PRODUCED IN GERMANY"
preserve 

local ele_type "windoffshore windonshore solar otherrenewable marine hydrowaterreservoir hydrorunofriverandpoundage hydropumpedstorage geothermal waste biomass nuclear fossilbrowncoallignite fossilcoalderivedgas fossilhardcoal fossilgas fossilpeat fossiloilshale fossiloil" 
egen load = rowtotal(`ele_type')

// Save dataset that contains updated "load" measure (based on total production)
keep TSO time_utc load
ren load load_supply
save load_entsoe_supply.dta, replace
restore


// Step length in line with technologies: consider 7 main technologies
egen step_0 = rowtotal(windoffshore windonshore solar otherrenewable marine hydrowaterreservoir hydrorunofriverandpoundage hydropumpedstorage geothermal)
gen step_1 = nuclear
egen step_2 = rowtotal(waste biomass)
gen step_3 = fossilbrowncoallignite
egen step_4 = rowtotal(fossilcoalderivedgas fossilhardcoal)
egen step_5 = rowtotal(fossilgas fossilpeat)
egen step_6 = rowtotal(fossiloilshale fossiloil)
mvencode step_*, mv(0) override // TSOs without these technologies


// Adjustment of oil plant in 50HZ as above
replace step_6 = 0 if TSO == "50Hertz"



*************************************************
// Figure 3 in the main text

// Marginal plants at 15 min interval
preserve
keep time_utc daily hourly omega_* marginal* TSO
drop if missing(time_utc)

tab marginal,gen(m_) sort

graph bar m_* , over(TSO) stack percentages ///
legend(label(1 "Natural gas") label (2 "Hard coal") label (3 "Nuclear") label (4 "Oil") label (5 "Brown coal") label (6 "Hydro") label (7 "Hydro - pumped") label (8 "Biomass"))

// graph export "../figures/marginal_plants.pdf", replace	

bysort TSO: tab marginal, sort 
restore




// DISTRIBUTION OF R - Q BY TSO
lab var omega_15 "{&omega} > 0,  15min"
hist omega_15 if omega_15 > 0, by(TSO) ///
	graphregion(lwidth(none) ilwidth(none)) plotregion(lwidth(none) ilwidth(none)) scheme(s1color)
// graph export "../figures/hist_omega_15_by_TSO.pdf", replace	

	
// COMPUTE EXPECTED VALUE OF DERIVATIVE OF AS WRT R
merge m:1 TSO time_utc using data_15min
keep if _merge == 3
drop _merge


// Need to merge the cluster for load 
ren daily date
ren hourly hour
merge m:1 TSO date hour using "data15min_cluster_TSO.dta", keepusing(kload)
keep if _merge ==3 
drop _merge


// DERIVATIVES OF REGRESSION MODEL WRT SOLAR (1_AS_regressions.do)
	gen dAS_dR = 5.174 +  (-0.00223) * 2 * solar_gen + ///
		( -9.83e-09) * 3 * solar_gen^2  + (0.000380) * load + ///
		 ( -6.76e-08) * load^2 + ( 0.000000189) * 2 * solar_gen * load if kload == 1 & TSO == "50Hertz"
		
	replace dAS_dR =  6.455 +  (-0.000391) * 2 * solar_gen + ///
		(-3.10e-08) * 3 * solar_gen^2  + ( 0.0000577) * load + ///
		 (-6.62e-08) * load^2 + (5.46e-08) * 2 * solar_gen * load if kload == 3 & TSO == "50Hertz"	
		
	replace dAS_dR =   18.44 +   (-0.00560) * 2 * solar_gen + ///
		(0.000000512) * 3 * solar_gen^2  + ( -0.00260) * load + ///
		  0.000000159 * load^2 + (0.000000185) * 2 * solar_gen * load if kload == 2 & TSO == "50Hertz"			  

	replace dAS_dR = 79.60 +  0.00101 * 2 * solar_gen + ///
		(0.000000101) * 3 * solar_gen^2  + (-0.00558) * load + ///
		 0.000000103 * load^2 + (-9.49e-08) * 2 * solar_gen * load if kload == 1 & TSO == "Amprion"
		
	replace dAS_dR =  -20.40 +  (-0.000250) * 2 * solar_gen + ///
		(0.000000442) * 3 * solar_gen^2  + (0.00184) * load + ///
		 -2.90e-08 * load^2 + (-0.000000117) * 2 * solar_gen * load if kload == 2 & TSO == "Amprion"	
		
	replace dAS_dR =  77.84 + 0.00388 * 2 * solar_gen + ///
		(-6.73e-09) * 3 * solar_gen^2  + (-0.00705) * load + ///
		  0.000000156 * load^2 + (-0.000000155) * 2 * solar_gen * load if kload == 3 & TSO == "Amprion"			  

	replace dAS_dR = 380.1 + (-0.000817) * 2 * solar_gen + ///
		(1.43e-08) * 3 * solar_gen^2  + (-0.0347) * load + ///
		  0.000000788 * load^2 + (3.51e-08) * 2 * solar_gen * load if kload == 1 & TSO == "TenneT"
		
	replace dAS_dR = -3.852 + 0.00105 * 2 * solar_gen + ///
		(-5.31e-09) * 3 * solar_gen^2  + (-0.0000909) * load + ///
		 1.98e-08 * load^2 + (-5.87e-08) * 2 * solar_gen * load if kload == 2 & TSO == "TenneT"	
		
	replace dAS_dR = 42.15 + 0.00161 * 2 * solar_gen + ///
		(-1.91e-08) * 3 * solar_gen^2  + (-0.00512) * load + ///
		  0.000000148 * load^2 + (-6.88e-08) * 2 * solar_gen * load if kload == 3 & TSO == "TenneT"			  

	replace dAS_dR = 130.8 +  0.00561 * 2 * solar_gen + ///
		(-1.23e-08) * 3 * solar_gen^2  + (-0.0329) * load + ///
		 0.00000202 * load^2 + (-0.000000616) * 2 * solar_gen * load if kload == 1 & TSO == "TransnetBW"
		
	replace dAS_dR = 40.25 + 0.0160 * 2 * solar_gen + ///
		(0.00000132) * 3 * solar_gen^2  + (-0.0217) * load + ///
		 0.00000274 * load^2 + (-0.00000375) * 2 * solar_gen * load if kload == 2 & TSO == "TransnetBW"	
		
	replace dAS_dR = 105.7 + 0.0130 * 2 * solar_gen + ///
		(-0.000000160) * 3 * solar_gen^2  + (-0.0314) * load + ///
		  0.00000223 * load^2 + (-0.00000150) * 2 * solar_gen * load if kload == 3 & TSO == "TransnetBW"			  

	  
	  
hist dAS_dR, by(TSO) xlabel(, grid)
// graph export "../figures/hist_dAS_dR_by_TSO.pdf", replace	


**********************************
// TABLE 1: MARGINAL BENEFITS

// LIMIT to solar_gen >0 for table with marginal benefits (standard deviations)
// keep if solar >0 
** keep all observations for reallocation exercise (robust)
**********************************

// COMPUTE EXPECTED VALUE OF MARGINAL COST
gen OC = marginal_cost * omega_15min

egen OC_weighted = total(OC), by(TSO)
sum OC_weighted if TSO == "50Hertz" 
sum OC_weighted if TSO == "Amprion" 
sum OC_weighted if TSO == "TenneT" 
sum OC_weighted if TSO == "TransnetBW" 
lab var marginal_cost "marginal cost ({c 0128}/MWh)"
hist marginal_cost, by(TSO) xlabel(, grid) color(blue%20) ///
	graphregion(lwidth(none) ilwidth(none)) plotregion(lwidth(none) ilwidth(none)) scheme(s1mono)
graph export "../figures/hist_OC_by_TSO.pdf", replace	


//SD for marginal cost
su marginal_cost if TSO == "50Hertz" 
su marginal_cost  if TSO == "Amprion" 
su marginal_cost  if TSO == "TenneT" 
su marginal_cost  if TSO == "TransnetBW" 


// COMPUTE EXPECTED VALUE OF MARGINAL EMISSIONS
gen avoid_emissions = marginal_emissions * omega_15min
egen avoid_emissions_weighted = total(avoid_emissions), by(TSO)
sum avoid_emissions_weighted if TSO == "50Hertz" 
sum avoid_emissions_weighted if TSO == "Amprion" 
sum avoid_emissions_weighted if TSO == "TenneT" 
sum avoid_emissions_weighted if TSO == "TransnetBW" 
hist marginal_emissions, by(TSO) xlabel(, grid)
// graph export "../figures/hist_avoid_emissions_by_TSO.pdf", replace	

drop if missing(time_utc)
// replace TSO = "50 Hertz" if TSO == "50Hertz"
saveold  "lambda_15min", replace


*******************************************************************************
local SCC 31.71   // IN EUR/tCO2


*******************************************************************************

tabstat avoid_emissions_weighted, by(TSO) s(mean N)

gen avoided_emissions_cost = `SCC' * avoid_emissions_weighted 
tabstat avoided_emissions_cost, by(TSO) s(mean N)

//SD
gen avoided_emissions_cost1 = `SCC' *marginal_emissions 
su avoided_emissions_cost1 if TSO == "50Hertz" 
su avoided_emissions_cost1  if TSO == "Amprion" 
su avoided_emissions_cost1  if TSO == "TenneT" 
su avoided_emissions_cost1  if TSO == "TransnetBW" 


tabstat dAS_dR, by(TSO) s(mean N)
*tabstat dAS_dR if kload==1, by(TSO) s(mean N)
*tabstat dAS_dR if kload==2, by(TSO) s(mean N)
*tabstat dAS_dR if kload==3, by(TSO) s(mean N)

su dAS_dR if TSO == "50Hertz" 
su dAS_dR if TSO == "Amprion" 
su dAS_dR if TSO == "TenneT" 
su dAS_dR if TSO == "TransnetBW" 


// OVERALL GAINS 
gen MB = OC_weighted + `SCC' * avoid_emissions_weighted - dAS_dR
tabstat MB, by(TSO) s(mean N)

// SD for total
gen MB1 = marginal_cost + avoided_emissions_cost1 - dAS_dR
su MB1 if TSO == "50Hertz" 
su MB1 if TSO == "Amprion" 
su MB1 if TSO == "TenneT" 
su MB1 if TSO == "TransnetBW" 
